Ageing and Relaxation in Glass Forming Systems 
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We propose that there exists a generic class of glass forming systems that have competing states 
(of crystalline order or not) which are locally close in energy to the ground state (which is typically 
unique). Upon cooling, such systems exhibit patches (or clusters) of these competing states which 
become locally stable in the sense of having a relatively high local shear modulus. It is in between 
these clusters where ageing, relaxation and plasticity under strain can take place. We demonstrate 
explicitly that relaxation events that lead to ageing occur where the local shear modulus is low 
(even negative), and result in an increase in the size of local patches of relative order. We examine 
the ageing events closely from two points of view. On the one hand we show that they are very 
localized in real space, taking place outside the patches of relative order, and from the other point 
of view we show that they represent transitions from one local minimum in the potential surface to 
another. This picture offers a direct relation between structure and dynamics, ascribing the slowing 
down in glass forming systems to the reduction in relative volume of the amorphous material which 
is liquid-like. While we agree with the well known Adam-Gibbs proposition that the slowing down 
is due to an entropic squeeze (a dramatic decrease in the number of available configurations), we do 
not agree with the Adam-Gibbs (or the Volger-Fulcher) formulae that predict an infinite relaxation 
time at a finite temperature. Rather, we propose that generically there should be no singular crisis 
at any finite temperature: the relaxation time and the associated correlation length (average cluster 
size) increase at most super-exponentially when the temperature is lowered. 



I. INTRODUCTION 



In this paper we attempt to provide generic answers 
to some long-standing questions regarding the spectacu- 
lar increase in relaxation times in typical structural glass 
formers as the temperature is reduced P, [H, H[ . In par- 
ticular we address the relation between structural and 
dynamical phenomena, the question of existence of a typ- 
ical length scale that is associated with the slowingdown 
[1, [H, [(j 0] , the nature of ageing in such systems j|[ and 
the issue of existence of a "true" glass transition temper- 
ature T g > where the relaxation time becomes infinite 
Q. We will base most of our comments on a careful 
analysis of one model system on which we performed ex- 
tensive simulations, but will provide supporting evidence 
also for a second model system and an additional ex- 
perimental system for which theory had been proposed 
recently. We will argue that these very different systems 
share some important characteristics, and we will risk a 
conjecture that these characteristics are generic for struc- 
tural glass formers. Most importantly, these systems 
share the characteristic of having, in addition to a crys- 
talline ground state, other states whose energy does not 
differ much from the ground state, at least for relatively 
small patches. The existence of such states leads first 
and foremost to the frustration of crystallization when 
the temperature is lowered. Glass formation is accom- 
panied by the creation of local patches, or clusters, of 
different nature, and these patches are locally stable in 
the precise sense of having high local shear moduli. To 
make this point crystal clear (pun not intended) we will 
introduce and study the notion of a local shear modulus 
for systems of this kind. The existence of the inhomo- 
geneity (or patchiness) is characterized by a typical scale, 
and we will demonstrate how this typical scale increases 



rapidly upon decreasing the temperature. Ageing and re- 
laxation events occur generically in the remainder of the 
system, in regions of low local shear moduli, as will be 
shown below. The squeezing out of the regions where re- 
laxations occur is the fundamental reason for the slowing 
down, and we connect the increase of typical scale with 
the increase of relaxation time. We will analyze individ- 
ual ageing events and show that they can be understood 
either as localized "excitation chains" [lj| in real space 
or transitions between local minima of the potential sur- 
face [ill • Finally we will argue that our picture docs not 
predict the existence of a 'true' glass transition, in the 
sense that infinite relaxation times are expected only as 
T — > (unless a very unlikely catastrophic crystallization 
event intervenes at a finite temperature). 

In Section [II] we present new results for the binary mix- 
ture glass-forming model, studying in a variety of ways, 
qualitative and quantitative, the spatial inhomogeneity 
that is crucial for glass formation. In Section IIIII we in- 
troduce the notion of "local shear modulus" , explain how 
to compute it in numerical simulations, and demonstrate 
the utility of this notion. We show that in the clusters of 
relative order the local shear modulus is typically high, 
whereas in between the clusters the local shear modulus 
can be low, even negative, signifying local mechanical in- 
stability. Ageing events take place spontaneously in the 
latter regions. In Sec. IIVI we consider the anatomy of 
the aging events, and show that they are very localized 
in real space. In terms of the potential energy landscape 
the aging events are mapped to transitions between lo- 
cal minima crossing saddles of order 1, (not necessarily 
always reducing the potential energy!). In Sect. IVl we 
discuss the generality of the picture emerging here with 
the help of other examples of glass formers. Sect lVII offcrs 
a summary and a discussion. 
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II. CLUSTERS AND HETEROGENEITY IN A 
BINARY GLASS FORMER 

The 2D system consists of an equimolar mixture of two 
types of point particles interacting via an inverse power 
potential [HEl: 

<Mr)=e(^)'\ a, b =1,2. (1) 

where r is the separation distance between particles of 
types a and b and n = 12. The particles have the same 
mass m, but half of the particles are 'small' with 'diam- 
eter' an = c, and half of the particles are 'large' with 
'diameter' 022 = 1.4c; the interaction between different 
kinds of particles is defined by C12 = (en + er 2 2)/2. We 
choose the units of mass, length and time as m, a and 
t = cy m/e , respectively. Other reduced units are the 
system energy U* = U/e, enthalpy H* = H/e, pressure 
P* = P<7 2 /e, reduced temperature T* = k B T/e, with 
Ub being Boltzmann's constant, and density p* = N/A*, 
where N is the particle number in the simulation box of 
dimensionless area A* = A /a 2 . 

We performed Molecular Dynamics simulations in the 
canonical (NVT) ensemble with N = 1024 particles in a 
square simulation box of side L* = s/A* with periodical 
boundary conditions. The equations of motion were inte- 
grated using a third order Gear predictor-corrector algo- 
rithm [TJ] with a time step St = 0.005-r. A constant tem- 
perature was preserved using a velocity rescaling method 
jl5l |. At each temperature the density was chosen in ac- 
cordance with the simulations results in an NPT ensem- 
ble as described in [l3[ with the pressure value fixed at 
P* = 13.5. 



A. Clusters and heterogeneity - qualitative picture 

1. Voronoi diagrams 

The Voronoi tessellation is an effective tool for the 
structural analysis of many-body systems. In the partic- 
ular case of two-dimensional systems the Voronoi poly- 
gon construction associates with every configuration of 
the particles a subdivision of position space into poly- 
gons, one per particle. The polygon associated with any 
particle contains all points closer to that particle than to 
any other particle. The edges of such a polygon are the 
perpendicular bisectors of the vectors joining the centers 
of neighboring particles. 

Any Voronoi construction obeys Euler's theorem : 
X = V — E + F, where x is the Euler characteristic, V, 
E and F are respectively the number of vertices, edges 
and faces. In an infinite two-dimensional system with 
periodic boundary conditions (on the torus X = 0), the 
average number of sides of the polygons is exactly 6 [l(| . 
This fixed average should be conserved under a local re- 
arrangement of the polygons. The elementary collective 




FIG. 1: A Tl process in the Voronoi tessellation (in red 
lines) . 

movement of particles which leads to local topological 
rearrangement of polygons, and known as Tl process, 
is shown in Fig. [TJ In a binary mixture there are two 
types of particles (small and large, or blue and red, re- 
spectively). Thus to achieve a mapping of the particle 
positions and the number of sides of the polygons in 
the Voronoi tessellation, we distinguish between polygons 
having a small or large particles in their center. Thus, a 
coloring scheme of cells will take into account not only 
the number of sides, but also the type of the particle in 
each cell. 



2. The competing phases 

At this point we would like to explain that the present 
system has one ground state, that of phase separated 
hexagons of small and large particles, but nearby there 
is another homogeneous state which is not too far in en- 
ergy. We expect that the ground state of the system 
under consideration should be homogeneous. The obvi- 
ous candidate structure is a phase-separated large and 
small particles (an example of Voronoi tessellation for 
subsystem of small or large particles is shown in the up- 
per panel of Fig. [J). But this is not the only possible 
homogeneous phase; It was shown 0, [l?], flq that the 
phase shown in the lower panel of Fig. [5] is possible for 
our binary systems in the range 1.18 < C22/C < 2. This 
phase consists of small particles in pentagons and large 
particles in heptagons; these are the polygons which in 
the theory of two-dimensional melting of one-component 
hexagonal crystals are referred to as 'defects' [l9| . This 
phase is characterized by pairs of small particles in pen- 
tagons lying in parallel lines [l7| ■ In one line all the pairs 
have the same orientation but can change it from line to 
line. Therefore, this phase is periodic in one direction 
and can be ordered or disordered in the other one. It is 
interesting to compare these two homogeneous states in 
terms of their energy and enthalpy, since we will argue 
below that they compete in the glassy state. 

The ground state of this model (at T = 0) is defined 
by minimum of the enthalpy H* = U* + P* ■ A* or the 
energy U* in an NPT or NVT ensemble respectively. In 
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lower temperatures dynamical relaxation slows down. A 
precise glass transition had not been identified in [l3j . In 
H, @ it was argued on the basis of statistical mechanics 
that there exists a typical length scale that grows super- 
cxponentially fast when the temperature decreases. As- 
sociated with this fastly growing scale there exists a fastly 
growing relaxation time, such that below a certain tem- 
perature the system is jammed for all practical purposes. 
Here we will shed further light on this phenomenon, re- 
lating it to the inhomogeneity seen with the bare eye in 
Fig.H 




FIG. 3: A Voronoi tesselation of a typical glassy state. Note 
the existence of clusters of large particles in hexagons, small 
particles in hexagons and patches of the phase shown in Fig. 
[2] Our point of view is that glass formers generically form 
clusters of such competing ground and close to ground states. 



FIG. 2: Upper panel: a Voronoi polygon construction for the 
hexagonal structure of one-component system. Lower panel: 
another phase of the binary mixture that can be stable at low 
temperatures, made from a Voronoi tessellation of heptagons 
with large particles and pentagons with small particles. 



III. CLUSTERS AND HETEROGENEITY- 
QUANTITATIVE PICTURE USING THE LOCAL 
SHEAR MODULUS 



the appendix we show that the ground state of our sys- 
tem is indeed the phase separated hexagons of small and 
large particles. The phase consisting of pentagons and 
heptagons is however close in energy, and certainly once 
formed will not be easy to deform to the phase separated 
ground state. 

The crucial comment is that in the super-cooled state 
this model exhibits clusters of all these phases. In Fig. 
[3] we show snapshots of the system at the temperature 
T* = 0.2. In the remainder of this paper we connect be- 
tween the heterogeneity that is caused by the existence of 
such clusters and the phenomenon of slow ageing in the 
glassy state. Ref. [l3[ found, using molecular dynamics 
simulations in the isothcrmal-isobaric (NPT) ensemble, 
that for temperature T > 0.5 the system is liquid and for 



In an NVT ensemble a given configuration can be char- 
acterized by components of the microscopic stress tensor 
which is defined by (see, e.g. pol]): 
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where pf is the a component of the dimensionless mo- 
mentum of particle i and rfj is the a component of the 
vector joining particles i and j. 

The diagonal components of this tensor correspond to 
density fluctuations and its trace determines the pres- 
sure, in accordance with the virial theorem: 
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FIG. 4: Trajectories from Molecular Dynamics simulations 
in the stress-energy plane, for the three temperatures T=0.2, 
T=0.3 and T=0.4 respectively. U* /N is the potential energy 
of the system per particle. Each point represents an average 
over a duration of 5r, connected with straight lines. The time 
of simulation is 10 4 r. Note the dramatic increase in the scale 
of the variations when the temperature increases. 



where d is the spatial dimension, [h] For interactions of 
the form (JTJ) the pressure is related to the total potential 
energy of the system (|A1[) : 



P* = p*T*+p 



dlv 



(4) 



The point is that an economic designation of a 'state' 
is given in the NVT ensemble by the pair of quantities 



energy and shear stress 

In Fig. H]we exhibit the actual trajectory of the Molec- 
ular dynamics simulations in the stress-energy plane for 
three different temperatures and for the same simulation 
time. We see that at T = 0.2 the trajectory hovers mostly 
around two distinct states with infrequent transitions be- 
tween them (only one transition in this run). This is 
a general observation for low temperatures: the system 
fluctuates around one "solid-like" state, but then jumps 
to another such "solid-like" state. The ease of transi- 
tions increases with temperature but also with the num- 
ber of particles. Here we have 1024 particles; in previous 
simulations using 256 particles [2l[ we did not observe 
any transition at T = 0.2. Changing the temperature 
to T = 0.3 but keeping the same simulation time one 
resolves seven-eight "solid-like" states. For T = 0.4 the 
trajectory now fills up an extended region in the stress- 
energy plane. As stated in [2l|, it appears that between 
'real' solid and 'real' liquid the system is locally not a 
"viscous fluid" . Rather, the trajectory concatenates rel- 
atively long period of time where the system behaves like 
a solid, interconnected by relatively short period of times 
where the system flows between these states. Clearly, 
a viscous fluid behaves very differently, responding to 
stress by a viscous flow, be it as sluggish as one wishes. 
Here, most of the time, the system does not respond to 
stress, except in the narrow corridors of flow which be- 
come rare when the temperature goes down and more 
common when the temperature warms up. Of course this 
does not mean that in the sense of long time averaging 
the notion of viscosity cannot be resurrected, but locally 
in time this is impossible. 



Since we have the system trajectory at our disposal, 
we can examine more closely the structural transforma- 
tion that is taking place when the system jumps from 
one "solid-like" state to the other. In Fig. O we show 
what happens by monitoring three characteristics of the 
system. One is simply the energy per particle as a func- 
tion of time cf. upper inset. We see that the energy 
fluctuates around a given value, and then jumps to a 
new state where the energy fluctuates at a lower value. 
A second characteristic is the average shear stress {cr xy ) 
(cf. lower inset) which exhibits a similar jump, except 
that the variance of its fluctuations before the transition 
appears considerably larger than the same variance after 
the jump. Most interestingly, the figure shows a window 
into the particle configuration itself, focusing on the local 
event that is responsible for the jumps in the insets. We 
see that the event is a concerted motion of six particles 
that change positions along a circular path; such events 
are referred to as "excitation chains" [lCj . At this point 
we want to explore what is the significance of this event, 
why it takes place where it docs, and what is it accom- 
plishing in terms of ageing and relaxation. To this aim 
we will study local mechanical properties. 
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FIG. 5: Details of an event of the type seen in Fig. pleading 
to a jump from one solid-like state to another. In the inset we 
show the energy per particle and the average shear stress as 
a function of time. Every point in these graphs is an average 
over 10 4 simulation steps. The figure itself is a window in the 
simulation box showing that the event is a circular motion 
of particles replacing each others position. The color code 
outside the event is red particles are large and black are small. 
In the circular event the red particle is large and all the rest 
are small 



A. Local shear moduli 

For inhomogcneous systems such as amorphous mate- 
rials it is important to study local properties in connec- 
tion with thermodynamical and structural characteris- 
tics. To measure local values we take our simulation box 
and subdivide it into smaller squares of length Ik = L* /2 k 
(2 < k < 4). We calculated the local stress and the local 
shear modulus following the definition of these quantities 
in [2^ . In the frame of this approach the local stress ten- 
sor ascribed to a subdomain to (1 < to < k 2 in the kth 
level) is given by 0, [H| : 



'a/3 



,T*5, 
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a/3 



E 



2 -It 



dr.; 



(5) 



where p* m is the particle number density in the subdo- 
main m, and is the length of the line segment of the 
vector joining particles i and j which is located inside a 
subdomain to (if this vector does not pass through this 
subdomain, (7™ = 0). The term q^/fij determines how 
different pairwisc interactions are apportioned to the lo- 
cal stress of a subdomain to, and includes contributions 
from pairs of segments located outside the subdomain. 
An average of <r^ over the entire area of the simulated 
system yields the usual bulk stress tensor ©. 

The local elastic modulus tensor is related to the in- 
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FIG. 6: Histograms of the local shear moduli for different 
temperatures (see the line-type in the inset) and for different 
partitions, from gross (k — 2, upper panel) to fine (k = 4, 
lower panel). Note the tendency to move the average to pos- 
itive values with lowering the temperature, with nevertheless 
a tail at T = 0.2 going into negative values of the local shear 
modulus. 



ternal stress fluctuations. The expression for the shear 
modulus per subdomain is defined by [22I [23j : 
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(«^ a p)~(a" a l )M), (6) 



where angular brackets mean averaging over configura- 
tions and the "infinite frequency shear modulus" (j,^ 
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[13, Ull P er subdomain is given by: 
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B. Results 

The results for the local shear modulus can be exhib- 
ited in a number of forms, each shedding light on another 
aspect of this useful measure. In Fig. [5] we show his- 
tograms of the local shear modulus, for different levels ol 
box resolution, from coarse to fine as one goes from the 
upper panel to the lowest. In general we see that as long 
as the temperature is high such that the system is fluid, 
the histograms peak around zero and are quite narrow. 
With T = 0.4 the histogram widens, but the mean still 
moves only a bit from zero. Not so for T = 0.2 which is 
deep in the glassy domain, where the histograms move 
to exhibit a positive average, but still showing wide vari- 
ations in the values of the local shear modulus. Note 
that with the increase of resolution we find more nega- 
tive values of the shear modulus, and these are indicative 
of sensitive regions in the system which are mechanically 
ready to exhibit ageing events. To clarify this point fur- 
ther we turn now to another way of presenting the same 
results, color coding real space to exhibit the high degree 
of spatial inhomogeneity of this system. 

Typical spatial results for the local shear modulus witl: 
spatial resolution k = 3 are shown in Fig. [Jj The loca 
shear modulus is color coded, showing positive and nega- 
tive values of the local shear modulus. Obviously, in the 
liquid (high temperatures) the distribution of local values 
averages to zero when summed over all the cells, whereas 
when the system is at low temperatures the average it 
positive, in agreement with the results of Ref. [2l[ . Also 
for higher temperatures the spread of values of the shean 
modulus (or the variance of its distribution) is smaller 
whearas the variance also increases upon decreasing the 
temperature. We can then expect that spontaneous age- 
ing events will take place where the local shear modulus 
is low. Indeed, for the system shown in Fig. [7J the firsl 
ageing event took place in the square (4,5), which is t 
region where the local shear modulus attains its lowesl 
value. We checked over many simulations that this is in- 
deed typical, and that looking at the distribution of local 
shear moduli one can guess very reliably where the next 
ageing event will take place. Of course, once the age- 
ing event occurred, the distribution of local shear moduli 
changes dramatically. In Fig. [5] we show the distribu- 
tion of local shear moduli after the event in square (4,5) 
took place. We see that the local shear moduli went up 
considerably in the region of the event; This is shown 
quantitatively in Fig. [9] where the difference between the 
local shear modulus in Figs. [5] and [7J is presented. We 
certainly need to understand what is precisely happening 
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FIG. 7: Color-coded local shear modulus before the ageing 
transition. Note that the local values can be either positive 
or negative, with positive values indicating local mechanical 
stability, whereas local negative values indicate mechanical 
instability. We expect relaxation (ageing) events to occur in 
the blue regions. 
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FIG. 8: Color-coded local shear modulus after the local re- 
laxation event occurred in box (4,5). 



here. A clue to that aim is obtained when we superpose 
on Fig. [7J the clusters of hexagonal patches of big par- 
ticles, see Fig. [TUJ We see that most of the changes in 
local shear modulus occur outside the clusters which are 
locally stable in the mechanical sense. 
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FIG. 11: ageing process. Shown is the potential energy per 
particle measured in the direct numerical simulation. Note 
the existence of five relaxation events. We analyze carefully 
the two events marked with ellipses. 



FIG. 9: Color-coded difference in the local shear modulus be- 
tween the two last figures, computed after the local relaxation 
event occurred in box (4,5). 
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FIG. 10: The same as Fig. [9] but with the clusters of large 
particles superposed to show that the main changes occur 
outside the clusters which are locally stable in the mechanical 
sense. 



IV. THE ANATOMY OF AGEING EVENTS 

In this section we examine the same localized events of 
ageing, or further relaxation steps towards equilibrium, 
but change the point of view from real-space to the po- 
tential energy landscape. We start with introducing the 
basic notions. 



A. Potential Energy Landscapes 

One popular way of discussing the dynamics of many 
particle systems is to visualize the state of the system as 
a point on the surface that the potential energy function 
U = f(x) draws in a Nd+ 1-dimensional space, where N 
is the number of particles and d the dimension of space 
(d = 2 in our examples) [2fj| . The evolution in time is de- 
scribed by the path of the coordinate vector x along the 
energy surface. Adding the thermal energy to the poten- 
tial energy elevates the point above the energy surface, 
but we will be interested in locating the paths associated 
with the relaxation events on the potential energy sur- 
face itself. At low temperatures the system is trapped 
in a minimum of this surface. In a crystallized system 
this minimum corresponds to the global minimum of the 
system - the ground state. In an amorphous system, 
however, this minimum is typically not a global one, but 
rather a local minimum which is separated from other 
minima of the system by moderate potential barriers. 
At temperatures greater than zero, we expect that, due 
to thermal fluctuations, the system will pass from one 
minimum to another via saddles. These saddles can be 
classified by their 'order', which by definition is the num- 
ber of negative eigenvalues of the Hessian matrix of U 
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at the saddle. At low tem per atures one encounters typi- 
cally saddles of order one [27| . Our aim here is to focus 
on spontaneous ageing events and map them onto trajec- 
tories on the energy surface. To do so we will take data 
from our numerical simulations, locate the ageing events, 
and for each point in time before and after the event find 
the nearest extrcmum of the potential. The methods 
used to find the extrema are standard, and are summa- 
rized for convenience in Appendix [Bj Unfortunately the 
numerics involved in this procedure is rather heavy, and 
for that reason in this section we will refer to simulations 
with 256 rather than 1024 particles and also somewhat 
higher temperatures where there is a larger number of 
ageing events, cf Fig. [U The insights achieved are of 
course independent of this numerical restriction. As an 
example consider Fig. [11] depicting the potential energy 
per particle as a function of time during an ageing pro- 
cess in a system of 256 particles. This figure contains 
evidence for five ageing events, and below we analyze in 
detail transitions occurring in the temporal proximity of 
the two events marked in red. These transitions are lo- 
calized, in the sense that the irreversible change in the 
configuration of the system involves only a few particles. 
In Figs. [T2l wc can see the actual changes in config- 
uration caused during these events. Here the particles 
positions are plotted before, during and after the event. 
The changed area consists of a mix of small and large 
particles. A temporary displacement occurs in a larger 
area around this region. This collective string-like pro- 
cess [28j allows for the permanent change to take place 
despite the high pressure. 

Fig. [13] shows the energy of the extrema found in the 
vicinity of these transitions. To find those we use a time- 
dependent trajectory from the molecular dynamics sim- 
ulation and for each configuration analyzed we seek the 
closest extrcmum of the potential energy surface. The 
vertical scale shows the potential energy of the system at 
the nearest extremum found using the LBFGS algorithm 
(see Appendix [B]). Most of the time the method locates 
a particular extremum with a specific energy, starting 
with a minimum, going through a saddle and ending at 
another minimum. Wc make sure that these extrema are 
indeed minima or saddles by finding the smallest eigen- 
value of the Hessian matrix at each extremum; it is pos- 
itive for the minima and unique negative in the saddle, 
meaning that the saddle is of order unity. 

Note that during the long residence near the first min- 
imum in Fig. [IS] there are instances when the system 
fluctuates to a point near a saddle; However the system 
remains trapped in the first minimum until the ageing 
transition takes place. 

The same method showed similar results for the first 
transition marked in Fig. 1111 i.e. a transition between 
two minima, with the transition state being a first order 
saddle point. It is therefore of interest to find the 're- 
action coordinate' for the transition, or the path on the 
potential surface which underlies the ageing process. It is 
very difficult to do so using the molecular dynamics sim- 




FIG. 12: (Color online) . Configurational changes during tran- 
sitions 1 and 2 in Fig. [11] respectively. Large symbols refer to 
large particles and small symbols to small particles. In circles 
we denote the positions of the particles before the transition, 
in x the positions at the saddle and in + the positions after 
the transition. 



ulation, and to achieve this we used the method of 'eigen- 
vector following' (see appendix for details) . Starting from 
the saddle, one takes a small step in one of the directions 
of the negative eigenvalue and then uses the eigenvector 
following algorithm to trace the path leading to the min- 
imum lying in this direction. This method works very 
well, resulting in the reaction coordinates shown in Fig. 
[T4l It is important to notice that the energies of the min- 
ima are indeed in agreement with the results of the pre- 
vious algorithm. We thus believe that the result shown 
connect indeed between the picture in time and the pic- 
ture in energy landscape, and the reaction coordinate is 
the one corresponding to the temporal events shown in 
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FIG. 13: Energy of nearest extrema as a function of time for 
the first and the second transitions in Fig. [TT] respectively. 
Not that in the upper panel we focus on an event in which 
the energy is increasing, and compare with Fig. 1141 



Fig. ED 



V. OTHER MODELS 

The emerging insight from the analysis presented in 
the pervious section is that the reason for the slowing 
down in glass forming systems is the generation of clus- 
ters of ground and close to ground states. These clusters 
have relatively high local shear moduli, since they are 
mechanically stable and can support locally a significant 
amount of shear without flowing. In between the clus- 
ters is where the ageing events take place, and indeed 
in those regions the local shear moduli are small, even 
negative, depending on the resolution of the calculation 
of the local shear modulus. When the temperature gets 
lower, the clusters increase in size, but there is no point 
of crisis where the relaxation time diverges unless T = 
or for some reason the system crystalizes. Note that the 
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20 40 60 80 100 120 140 160 180 200 220 
Algorithm Steps 



FIG. 14: "Reaction coordinate" for the two marked transi- 
tions in Fig. 1111 These graphs are obtained by starting at 
the saddle point and moving in the direction of the (unique) 
unstable eigenvector. Note that in general saddles are not 
guaranteed to be of order 1, and may have more than one 
unstable direction. 

lack of such finite temperature crisis was shown rigor- 
ously for this particular model of binary mixture in [t| . 
One question that we need to ask is whether this behav- 
ior is generic to glass formers or is it a peculiarity of the 
present model of binary mixture. 

In the context of our own work we have considered 
recently two other, very different, models of glass formers, 
i.e. the Shintani-Tanaka model [29[ and dry glycerol [3fj| . 
In both systems it was demonstrated, using quantitative 
theory, that the scenario discussed here appears general, 
independent of the very different details characterizing 
these systems @, Hl[ • Thus for example we could show 
in the context of the Shintatni-Tanaka model that there 
exists a typical sale £ that dominates the mean relaxation 
time as a function of the temperature. In other words, the 
T a relaxation time in this model could be quantitatively 
determined by a formula reading 

r a = t exp(A£/T) , (8) 

where A is a temperature-independent dimensional con- 
stant and To is the microscopic attempt time. A similar 
result was demonstrated for the present binary model in 
[j| . The non-Arrhenius behavior of this formula is due to 
the strong increase of £ when the temperature decreases, 
due to the increase in the size of the clusters, but without 
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any singularity at any temperature T > 0. For the glyc- 
erol one could go even further, predicting not only the 
mean relaxation time but also computing the functional 
form of the dielectric response function in good relation 
to experiments. Here one argued that each cluster of s 
molecules contributes an exponential decay of its dipole 
moment with an ,s-dcpendcnt decay constant. Once aver- 
aged over all the clusters, the resulting response function 
is strongly non-exponential as is indeed observed. The 
main point is that the long relaxation times are again 
due to the large clusters (smaller cluster may contribute 
(3 peaks when the conditions are right, cf. [31|). As long 
as the system does not crystalize there is no mechanism 
for a crisis of the type predicted by the Vogel-Fulcher 
formula or the Kauzmann argument. For more detailed 
discussion of this important issue the reader is referred 
to 1. 



VI. SUMMARY AND DISCUSSION 

In summary, we focused on generic ageing events in 
our glassy system, to understand the relation between 
slowing down and spatial inhomogeneity. We have shown 
that the inhomogeneity can be seen in two complimentary 
ways, one by observing the clusters of competing ground 
and close to ground states, and the other by measuring 
the local shear modulus. The latter quantity fluctuates 
from place to place, tending to be high and positive in the 
presence of clusters and low and even negative outside the 
clusters. Ageing events occur in these mechanically un- 
stable regions, and they result in a local re-organization, 
usually increasing the area spanned by clusters. We ex- 
amined closely the nature of these ageing events, and 
showed again that they can be described in two comple- 
mentary forms. One is a localized chain of movements, a 
so-called "excitation chain" (accompanied by a larger re- 
versible collective motion) and the second is a transition 
between two inherent states, or local minima in the po- 
tential surface of the system, with the transition crossing 
a saddle point of order unity. Note that we could dif- 
ferentiate between the very localized permanent change 
and the "string-like" event that is in part reversible. 

We proposed, on the basis of the present model of a 
binary mixture and other systems that were analyzed 
recently in the same spirit, that the results obtained here 
are generic, and that there exists a wide class of glass 
forming systems whose theoretical understanding calls 
for an analysis of its distribution of clusters of competing 
ground and close to ground states; in a future publication 
we will show that the temporal properties of such systems 
can be completely understood in these terms. 



APPENDIX A: THE GROUND STATE 

The average potential energy per particle in a binary 
mixture with the interactions defined by the potential ((T|) 



is given by: 



N 2 ■ N 2—> I n . 



12 



(Al) 



where is the distance between particles i and j. 

For the equimolar system of N — > oo separated small 
and large particles (see the upper panel of Fig. [J) the in- 
teractions between the two species is neglible (o(l/VN)) 
and the total energy (|A1[) can be written as: 



U* 
~N 
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(A2) 



where ru and ri2 are the distances between small and 
large particles respectively. 

The area of a hexagon is -^r 2 , where r is the distance 
between particles. The dimensionless total area A* = 
l/p* per particle of the system of separated small and 
large particles is given by: 



A* 
~N 



73 
4 



(A3) 



The enthalpy depends on two variables ru and r22, so 
the minimum of this function is defined by two equations 
dH* /dm — and dH* /dr 2 2 = 0. Solutions of these 
equations and evaluated values for the energy, enthalpy 
and density are displayed in Tab. |TJ 

TABLE I: a-NPT ensemble, b-NVT ensemble. 



Phase 


Sep. a 


Mixt. a 


Sep. b 


Mixt. b 


ru 


1.031 


1.003 


1.039 


0.995 


r22 


1.376 


1.473 


1.386 


1.462 


ri2 




1.187 




1.179 


p * 


13.5 


13.5 


12.2 


14.9 


P* 


0.781 


0.77 


0.77 


0.781 


U*/N 


2.881 


2.922 


2.648 


3.179 


H*/N 


20.168 


20.453 


18.752 


22.257 



In the configuration where small and large particles 
are mixed (the lower panel of Fig. [2]) each small particle 
interacts with one small particle and four large particles. 
In turn, each large particle interacts with three neighbor- 
ing large particles and four neighboring small particles. 
Therefore, the general expression (|A1|) is simplified to: 



U* 
~N 



on 

m 



012 

ri2 



Q"22 
7"22 



(A4) 



where in addition to (|A2|) the contribution of interactions 
between small and large particles separated by r\2 are 
taken into account. 

Any configuration of homogeneously mixed small and 
large particles can be constructed using 'boat' unit [l8[ 
(see the lower panel of Fig. ^ which consists of one 'thin' 
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and three 'fat' rhombi. There arc 1/4JV 'thin' rhombi and 
3/4iV 'fat' rhombi in the system, so the area per particle 
is given by: 



^ = l-(A^ + 3-A* F ), 



(A5) 



where is the area of a thin rhombus and A* F that of 
a fat one. 

It follows from the 'boat' design that the distances be- 
tween particles are not independent. It is useful to intro- 
duce a new variable a, the small angle of the 'fat' rhomb, 
so that the large angle of the 'thin' rhomb is 2tt — 3a. All 
rhombi have sides of the same size r\2 and the distances 
between small-small and large-large particles can be ex- 
pressed as: 



rn = -2 • r 12 ■ cos(-a) 
,1 . 

r 22 = 2 • ri2 • sin(-a) 

The area per particle of the mixed system 
function of these variables is given by: 

A* 1 

— = - • r\ 2 ■ (-sm(3a) + 3 • sin(a)) 
= r\ 2 ■ sin 3 (a). 



(A6) 
3) as a 

(A7) 



Substitution of (|A6p into (|A4[) yields the following ex- 
pression for the energy: 
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The minimum of the enthalpy in this case is defined 
by equations dH*/dri 2 = and dH* /da = 0. The so- 
lution for the angle a depends only on the ratio (T22/C11 
and for the system under consideration a = 76.65°. The 
interparticle distances, energy, enthalpy and density are 
shown in Table [T] 

In the NVT ensemble it is necessary to find the min- 
imum of the energy given by (|A2[) or (|A8|) under the 
constraints (|A3|) or (|A7|) . respectively. The calculation 
results are also presented in Tab|TJ 



APPENDIX B: LOCATING EXTREMA 

Finding a local minimum of a function of several vari- 
ables is a common problem in applied math. There are 
several methods for doing this, some of them like the 
Newton-Raphson method are quite simple. Finding a 
saddle point, on the other hand, is a more subtle prob- 
lem since it is not clear if the appropriate direction to go 
is up or downhill. In our work we have used two different 
methods to overcome this problem: 



1. Square gradient minimization 

One way of finding saddle points is by finding a local 
minimum of the squared modulus of the gradient of the 
original potential function: 

N 2 , „ s 2 

w(x) = \vu(x)\>=j2Y l l£ L ) • ( B1 ) 

i=l a=X \ % < a ' 

Since W(x) is nonnegative, at the absolute minima 
W(x) = 0, which implies Vf7 = (a saddle point). An 
efficient way to minimize such functions is the Newton- 
Raphson method. At each algorithm step = x n + 
5x n the following equation should be satisfied: 



VU(x n + 5x n ) = . 
This equation is solved to first order: 

VU(x n ) + 8x n V 2 U{x n ) = . 



(B2) 



(B3) 



Inverting for 8x n one finds 

Sx n = ~X7 2 U{x n )VU{x n ) = -HVU{x n ) . (B4) 

We have used a publicl y av ailable numerical algorithm 
- the LBFGS algorithm ( |32l. l33l]) which uses an approx- 
imate version of the Newton method. The biggest draw- 
back in minimizing the function W(jx.) instead of mini- 
mizing U (x) directly is that the minimization often con- 
verges to a local minimum of W(x) , which is not a saddle 
point of J7(x). Threfore, when minimizing a molecular 
dynamics realization of a system travelling in the energy 
landscape, not all of the instantaneous states converge to 
a minimum or a saddle of V(x) and some information is 
lost. 



2. Eigenvector following 

In order to demonstrate the topological connectedness 
of two minima separated by a saddle we have used the 
Eigenvector Following Method following the algorithm 
as described by Wales and coworkers [H, [H, Ull (the 
original idea for the algorithm was proposed by Cerjan 
and Miller [37|). At each iteration a step An; is proposed, 
which in the base that diagonalizes the Hessian (at the 
specific point in the algorithm path) is [H, [H| 



Aa; u — S u 



2£p 



(B5) 



where are the eigenvalues of the Hessian and <? M are 
the components of the gradient in the diagonal base (Ax p 
is set to for the directions where L = 0; i.e. uniform 



displacements). The sign 
below. Note that as — * 0, 



±1 is chosen as explained 



Ax,. = - 



h.. 



0(gl), 
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where the first term is the Newton-Raphson step. 

If = — 1 the algrithm converges to a minimum along 
jj,. Therefore, starting from a saddle of order one, setting 



all of the Sft to —1 we were able to reach the 2 minima 



separated by this saddle. 
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